### Functions for calculating cross elas matrix ###
# based upon counterpart functions in CF_MMX_functions.R BUT there are some differences

###ownership matrix n.models X n.models
build.co <- function(firm) {
  co <- matrix(0,nrow=length(firm),ncol=length(firm))
  for(i in 1:length(firm)) {
    for(j in 1:length(firm))
      if(firm[i]==firm[j] & i!=j) co[i,j] = 1  
  }
  return(co)
}

#Data generating process: based on DGP.exog() in CF_MMX_functions.R
DGP.BLP <- function(year.chosen,vdgp){
  #First adapt the heterogeneity to the chosen v
  if(vdgp==2){
    y.sd <- 0 # to get rid of price het. 
    alpha <- -BLP.meanownelas*ybar/BLP.priceadjust # this is adjusted to mimic the average own elasticity average
              
    }
  if(vdgp==1){
    y.sd <- 0 # to get rid of price het. 
    alpha <- -BLP.meanownelas*ybar/BLP.priceadjust # this is adjusted to mimic the average own elasticity average
    b.sds <- c(0,0,0,0,0)
    names(b.sds) <- param[param_name %like% "sigma"]$xvar
  }
  #
  cat("Inside DGP.BLP we have alpha = ",alpha,"y.sd = ",y.sd,"b.sds = ",b.sds, "\n")
  ##Demand-side
  ##income parameters
  grid.hh  <- data.table(hh=1:n.hh,year=year.chosen)
  mu.ly <- y.mn[year==year.chosen]$V2
  sigma.ly <- y.sd
  if(sequences %in% c("halton","sobol")) {
    if(sequences=="halton") {
      y.h.std <- randtoolbox::halton(n.hh,dim=1,normal=TRUE,init=TRUE) 
    } else y.h.std <- randtoolbox::sobol(n.hh,dim=1,normal=TRUE,scrambling = TRUE, seed = seedn) 
    grid.hh[, y.h := exp(mu.ly +sigma.ly*y.h.std) ] 
    grid.hh[, alpha.h := alpha/y.h]
    if(sequences=="halton") {
      b.h.std <- randtoolbox::halton(n.hh,dim=5,normal=TRUE,init=FALSE) 
    } else b.h.std <- randtoolbox::sobol(n.hh,dim=5,normal=TRUE,scrambling = TRUE) 
    grid.hh[, b.h.const :=  b.means[["const"]]+b.sds[["const"]]*b.h.std[,1]]
    grid.hh[, b.h.hpwt :=  b.means[["hpwt"]]+b.sds[["hpwt"]]*b.h.std[,2]]
    grid.hh[, b.h.air :=  b.means[["air"]]+b.sds[["air"]]*b.h.std[,3]]
    grid.hh[, b.h.mpd :=  b.means[["mpd"]]+b.sds[["mpd"]]*b.h.std[,4]]
    grid.hh[, b.h.space :=  b.means[["space"]]+b.sds[["space"]]*b.h.std[,5]]
  } else {  # code using rnorm draws
  grid.hh[, y.h := exp(rnorm(n=n.hh,mean=mu.ly,sd=sigma.ly)) ] 
  grid.hh[, alpha.h := alpha/y.h]
  grid.hh[, b.h.const :=  rnorm(n=n.hh,mean=b.means[["const"]],sd=b.sds[["const"]])]
  grid.hh[, b.h.hpwt :=  rnorm(n=n.hh,mean=b.means[["hpwt"]],sd=b.sds[["hpwt"]])]
  grid.hh[, b.h.air :=  rnorm(n=n.hh,mean=b.means[["air"]],sd=b.sds[["air"]])]
  grid.hh[, b.h.mpd :=  rnorm(n=n.hh,mean=b.means[["mpd"]],sd=b.sds[["mpd"]])]
  grid.hh[, b.h.space :=  rnorm(n=n.hh,mean=b.means[["space"]],sd=b.sds[["space"]])]
  }
  demandside <- data.frame(grid.hh[,list(year=year,alpha.h,b.h.const,b.h.hpwt,b.h.air,b.h.mpd,b.h.space,y.h)])  
  #
  ##Supply-side
  supplyside <- BLP[year == year.chosen]
  supplyside[,xi_m := 1] # starting value
  supplyside[,alpha.base := alpha] # in order to compute expected ownelas later
  #Mahalanobis distance caculations:
  if(y.sd==0) MW <- BLP[year == year.chosen,.(m.text=m,hpwt,air,space,mpd)] else {
    MW <- BLP[year == year.chosen,.(m.text=m,hpwt,air,space,mpd,p)]  # keep p when alpha het
  }
  #
  MW[,m := 1:.N]
  ML <- melt(MW,id=c("m","m.text")) # ML is the long version of MW (models, wide)
  SD <- ML[,.(S_v=sd(value)),by=.(variable)]
  X <- CJ(m=unique(MW$m),j=unique(MW$m),variable =SD$variable)  
  X <- merge(X,ML,by=c("m","variable"))
  setnames(X,c("value","m.text"),c("X_vm","m.text_m"))
  X <- merge(X,ML,by.x = c("j","variable"),by.y=c("m","variable"))
  setnames(X,c("value","m.text"),c("X_vj","m.text_j"))
  X <- merge(X,SD,by="variable")
  X[, sumterm := (X_vm-X_vj)^2/(S_v^2)]
  X[S_v==0, sumterm := 0] # AC has no variation in 1971
  MD <- X[,.(maha_dist = sqrt(sum(sumterm))),by=.(m,m.text_m,j,m.text_j)]
  #
  return(list(models=supplyside,households=demandside,distance=MD)) #comprises three data frames, 
  #one w/ n.models  rows and other w/ n.hh X n.countries rows, and final one with mahalanobis dist
}

#
# Computing ownelas ====
#
xi.BLP.mlog <- function(dat,verbose=FALSE) { 
  #set verbose = TRUE if you want to see the Iter and Residual  for the FPITER.
  source("Rfiles/eqbm_MLOG_subfunctions.R",local=TRUE)
  # xi update formula
  xi.update <- function(xi.loop,nu=0.3) {
    bxap <- b.h %*% t(x)-alpha.h %*% t(p)
    delta.old <- bxap + rep(1,n.hh) %*% t(xi.loop)
    s.predicted <- colMeans( exp(delta.old)/(1+rowSums(exp(delta.old))) )
    delta.new <- delta.old  + rep(1,n.hh) %*%t(log(s.true))  - rep(1,n.hh) %*%t(log(s.predicted))
    xi.new <- colMeans(delta.new - bxap)
    return(as.vector( nu*xi.new + (1-nu)*xi.loop) )
  }
  #
  # start of elasticities computation #  
  df.m <- dat$models
  #new
  xvars <-  names(b.means)
  x <- as.matrix(df.m[,..xvars])  # df.m is actually a data.table
  xi.init <- df.m$xi 
  p.true <- df.m$p
  s.true <- df.m$s
  n.models <- length(s.true)
  #
  df <- dat$households
  #new
  bvars <- paste0("b.h.",xvars)
  b.h <- as.matrix(df[,bvars])  # df is actually a data.frame
  #
  alpha.h <- df$alpha.h
  y.h <- df$y.h
  Y <- sum(df$y.h)
  #FPI 1: to get the xi
  co.own <- co.own.mat
  p <- p.true
  FPIresults <- fpiter(xi.init,fixptfn=xi.update,nu=1,control=list(trace=verbose,maxiter=MAXIT))
  n.loops.xi  <- FPIresults$fpevals
  xi <- FPIresults$par # model level price vector
  s <- colMeans(prob.mat(p))
  cat(paste0("Xi converged: ",n.loops.xi<MAXIT,"\n"))
  cat(paste0("Market shares match for year ",df.m[1]$year,": ", all.equal(s,s.true),"\n"))
  #collect model level variables
    m.eqvar <- data.frame(df.m$m ,xi.predicted_m=xi) # model level output for model m
    return(m.eqvar)
    }

elas_price <- function(p.change,dat,xi.dat,m.chosen){
  source("Rfiles/eqbm_MLOG_subfunctions.R",local=TRUE)
  # start of elasticities computation #  
  df.m <- dat$models
  #new
  xvars <-  names(b.means)
  x <- as.matrix(df.m[,..xvars])  # df.m is actually a data.table
  xi <- xi.dat$xi.predicted 
  p.true <- df.m$p
  s.true <- df.m$s
  n.models <- length(s.true)
  #
  df <- dat$households
  #new
  bvars <- paste0("b.h.",xvars)
  b.h <- as.matrix(df[,bvars])  # df is actually a data.frame
  #
  alpha.h <- df$alpha.h
  y.h <- df$y.h
  Y <- sum(df$y.h)
  #
  #
  p <- df.m$p
  m <- df.m$m
  p[m==m.chosen] <- p[m==m.chosen]*(1 + p.change) 
  s <- colMeans(prob.mat(p))
  pelas.m <- ownelas(prm=prob.mat(p),share=s,price=p) # model level own price elasticities
  superelas.m <- superelas(prm=prob.mat(p),share=s,price=p) # model level own super elasticity
  return(data.frame(p.change=p.change,m=m.chosen,p=p[m==m.chosen],eps=-pelas.m[m==m.chosen],E=superelas.m[m==m.chosen],s=s[m==m.chosen]))
}



superelas_mean <- function(dat,xi.dat){
  source("Rfiles/eqbm_MLOG_subfunctions.R",local=TRUE)
  # start of elasticities computation #  
  df.m <- dat$models
  #new
  xvars <-  names(b.means)
  x <- as.matrix(df.m[,..xvars])  # df.m is actually a data.table
  xi <- xi.dat$xi.predicted 
  p.true <- df.m$p
  s.true <- df.m$s
  n.models <- length(s.true)
  #
  df <- dat$households
  #new
  bvars <- paste0("b.h.",xvars)
  b.h <- as.matrix(df[,bvars])  # df is actually a data.frame
  #
  alpha.h <- df$alpha.h
  y.h <- df$y.h
  Y <- sum(df$y.h)
  #
  p <- df.m$p
  m <- df.m$m
  s <- colMeans(prob.mat(p))
  pelas.m <- ownelas(prm=prob.mat(p),share=s,price=p) # model level own price elasticities
  superelas.m <- superelas(prm=prob.mat(p),share=s,price=p) # model level own super elasticity
  pte.sp.m <- (pelas.m + 1)/(pelas.m + 1 - superelas.m)
  #take mean over all models (alternatively could be just the foreign models)
  return(data.table(eps.mean = mean(pelas.m),E.mean = mean(superelas.m),pte.mean = mean(pte.sp.m)))
}





#MLOG version of richsubs, now  includes an argument for whether we 
#intend to do only richsubs (fast and no convergence issue), or CF
#which means solving for mc, and verifying that we get back the true prices.
#
rich.subs.BLP.mlog <- function(dat,verbose=FALSE,objective="richsubs") { 
  #set verbose = TRUE if you want to see the Iter and Residual  for the FPITER.
  source("Rfiles/eqbm_MLOG_subfunctions.R",local=TRUE)
  # xi update formula
  xi.update <- function(xi.loop,nu=0.3) {
    bxap <- b.h %*% t(x)-alpha.h %*% t(p)
    delta.old <- bxap + rep(1,n.hh) %*% t(xi.loop)
    s.predicted <- colMeans( exp(delta.old)/(1+rowSums(exp(delta.old))) )
    delta.new <- delta.old  + rep(1,n.hh) %*%t(log(s.true))  - rep(1,n.hh) %*%t(log(s.predicted))
    xi.new <- colMeans(delta.new - bxap)
    return(as.vector( nu*xi.new + (1-nu)*xi.loop) )
  }
  #
  # (implied) cost update formula
  cost.update <- function(cost,nu) {
    #
    cannib2 <- function(prm,price)  {
      D.cross <- build.DeltaX(prm)
      return(D.cross %*% (price-cost))  # cost to be taken from function's env
    }
    #market shares
    pim <- prob.mat(p) # price vector in the environment.
    # apply the first order condition for profit-max. # c -(s+r)/(-d)
    mc.star <- p + (colMeans(pim)+cannib2(pim,p))/ownderiv(pim)  # cannib needs to take cost as an argument
    # weight the nu price
    return( as.vector( nu*mc.star +(1-nu)*cost ))
  }
  #  
  # n.hh length vector of sales-shares for model m. (nothing to do with ownership matrix)
  omega <- function(prm,m) {
    (prm[,m])/sum(prm[,m])  
  }
  #  
  #household variance in market share of m
  wtd.var.prm <- function(prm) {
    wt <- rep(1/n.hh,n.hh)
    mu <- colSums(prm*wt)  
    colSums(wt*(t(t(prm)-mu))^2) 
  }
  # crosselas : based on build.DeltaX, a matrix of cross price elasticities, changing price of m (row), impact on j(cols)
  crosselas.mat <- function(prm,price) {
    # should be simplified using this pim <- prob.mat(price)
    Delta <- matrix(0,nrow=n.models,ncol=n.models)
    for(m in 1:n.models) {
      for(j in 1:n.models) {
        Delta[m,j] = sum(alpha.h*prm[,m]*omega(prm,j))  # epsilon_{jm} * Omega_{jm} in theory
      }
    }
    return(price*Delta) # multiplies the matrix m rows j columns, by prices of m
  }
  #
  cov.mat <- function(prm,price) {
    # should be simplified using this pim <- prob.mat(price)
    Delta <- matrix(0,nrow=n.models,ncol=n.models)
    for(m in 1:n.models) {
      for(j in 1:n.models) {
        Delta[m,j] = (n.hh*sum(prm[,m]*prm[,j]))/(sum(prm[,j])*sum(prm[,m])) - 1    #  covariance divided by s_j*s_m
      }
    }
    return(Delta) 
  }
  #
  # start of elasticities computation #  
  df.m <- dat$models
  #new
  xvars <-  names(b.means)
  x <- as.matrix(df.m[,..xvars])  # df.m is actually a data.table
  xi.init <- df.m$xi 
  p.true <- df.m$p
  s.true <- df.m$s
  n.models <- length(s.true)
  #
  df <- dat$households
  #new
  bvars <- paste0("b.h.",xvars)
  b.h <- as.matrix(df[,bvars])  # df is actually a data.frame
  #
  alpha.h <- df$alpha.h
  y.h <- df$y.h
  Y <- sum(df$y.h)
  #FPI 1: to get the xi
  co.own <- co.own.mat
  p <- p.true
  FPIresults <- fpiter(xi.init,fixptfn=xi.update,nu=1,control=list(trace=verbose,maxiter=MAXIT))
  n.loops.xi  <- FPIresults$fpevals
  xi <- FPIresults$par # model level price vector
  s <- colMeans(prob.mat(p))
  cat(paste0("Xi converged: ",n.loops.xi<MAXIT,"\n"))
  cat(paste0("Market shares match for year ",df.m[1]$year,": ", all.equal(s,s.true),"\n"))
  #IF OBJECTIVE IS TO CALIBRATE CF:
  if(objective!="richsubs"){
  # FPI 2: infer the marginal costs from prices and the demand curvature (a la BLP)
  FPIresults <- fpiter(p,fixptfn=cost.update,nu=1,control=list(trace=verbose,maxiter=MAXIT))
  n.loops.mc  <- FPIresults$fpevals
  mc <- FPIresults$par
  cat(paste0("Marginal costs converged: ",n.loops.mc<MAXIT,"\n"))
  #
  # FPI 3: verify we get back the correct prices.
  p.monop.comp <- mc + 1/mean(alpha.h)
  FPIresults <- fpiter(p.monop.comp,fixptfn=price.update,nu=.3,control=list(trace=verbose,maxiter=MAXIT))
  n.loops.p  <- FPIresults$fpevals
  p <- FPIresults$par
  cat(paste0("Prices converged: ",n.loops.p<MAXIT,"\n"))
  cat(paste0("Prices match for year ",df.m[1]$year,": ", all.equal(p,p.true),"\n"))
  #
  }
  #those two are common to both objectives
  pelas.m <- ownelas(prm=prob.mat(p),share=s,price=p) # model level own price elasticities
  superelas.m <- superelas(prm=prob.mat(p),share=s,price=p) # model level own super elasticity
  ##also create averages and standard deviations over households
  #an n.hh length vector of average elasticities
  lownelas.hh <- log(-hhelas(p))
  lownelas.mean.hh <- mean(lownelas.hh)
  lownelas.sd.hh <- sd(lownelas.hh) #
  wvp <-  wtd.var.prm(prob.mat(p)) # model level market share vector
  #collect model level variables
  if(objective!="richsubs"){
  m.eqvar <- data.frame(df.m ,xi.predicted_m=xi,s.predicted_m=s,wvp_m=wvp,own_elas_m=pelas.m,super_elas_m =superelas.m,mc_m=mc,p.predicted_m=p) # model level output for model m
  } else {
  m.eqvar <- data.frame(df.m ,xi.predicted_m=xi,s.predicted_m=s,wvp_m=wvp,own_elas_m=pelas.m,super_elas_m =superelas.m) # model level output for model m
  }
  setnames(m.eqvar,old=c("m","air","hpwt","mpd","space","p","s")
           ,new=c("m.text","air_m","hpwt_m","mpd_m","space_m","p_m","s_m"))
  m.eqvar <- data.frame(m.eqvar, m = as.numeric(row.names(m.eqvar)))
  #
  cross.elas <- as.data.table(crosselas.mat(prm=prob.mat(p),price=p)) # model-level cross price elasticities matrix
  cov.prob <- as.data.table(cov.mat(prm=prob.mat(p),price=p)) # model-level matrix of covariances divided by s_j 
  #
  #reshape long
  cross.elas[, m := as.numeric(row.names(cross.elas))]
  cross.elas.long <- melt(cross.elas,measure.vars=patterns("^V"),value.name="cross_elas",id="m")
  cross.elas.long[,j := as.numeric(sub("V","",variable))]
  cross.elas.long[,variable := NULL]
  #
  cov.prob[, m := as.numeric(row.names(cov.prob))]
  cov.prob.long <- melt(cov.prob,measure.vars=patterns("^V"),value.name="cov_prob",id="m")
  cov.prob.long[,j := as.numeric(sub("V","",variable))]
  cov.prob.long[,variable := NULL]
  #
  #merge ownelas
  cross.elas.xdiff.ownelas <- merge(cross.elas.long,m.eqvar,by="m",all.x=TRUE)
  #add the scaled covariances
  cross.elas.xdiff.ownelas <- merge(cross.elas.xdiff.ownelas,cov.prob.long,by=c("m","j"),all.x=TRUE)
  #add the MDs
  cross.elas.xdiff.ownelas <- merge(cross.elas.xdiff.ownelas,dat$distance,by=c("m","j"),all.x=TRUE)
  #add the n.loops
  if(objective!="richsubs"){
    n.loops <- cbind(n.loops.xi,n.loops.mc,n.loops.p)
  } else {
    n.loops <- n.loops.xi
  }
  #
 # return(cbind(cross.elas.xdiff.ownelas,n.loops))
  return(cbind(cross.elas.xdiff.ownelas,n.loops,lownelas.mean.hh,lownelas.sd.hh))
  }

#Monte function
monte.func.rich.subs.mlog <- function(repi,year.monte,objective="richsubs",vdgp) { # monte carlo stage
  data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=vdgp) #gen exog data
  my.elas <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective=objective) 
}

#Monte rich subs moments
monte.func.rich.subs.mlog.moments <- function(repi,year.monte,vdgp) { # monte carlo stage
  data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp) #gen exog data
  DT <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective="richsubs") 
  res_loglin <- feols(log(cross_elas)~ maha_dist | m+j,DT[m !=j])  
  crosselasmaha <- res_loglin$coefficients[[1]] 
  meanmaha <- DT[m!=j, mean(maha_dist)] # IQR = 2 for 1990.
  IQR <- DT[m!=j, IQR(maha_dist)] # IQR = 2 for 1990.
  covprob <- DT[m!=j, mean(cov_prob)] # 
  exvar <- DT[m==j, mean(wvp_m/(s_m*(1-s_m)))] # excess variance, mn = 0.015
  DTcross <- DT[m!=j]
  DTcross <- DTcross[,.(sd.cross_elas_m=sd(cross_elas),mn.cross_elas_m=mean(cross_elas)), by=.(m)]
  mnsd <- DTcross[,mean(sd.cross_elas_m)] # mean of s.d. of cross elasticities = 0.021
  mnmn <- DTcross[,mean(mn.cross_elas_m)] 
  mncv <- DTcross[,mean(sd.cross_elas_m/mn.cross_elas_m)] 
  return(data.table(repi=repi,crosselasmaha=crosselasmaha,meanmaha=meanmaha,IQR=IQR,exvar=exvar,covprob=covprob,mnsd=mnsd,mnmn=mnmn,mncv=mncv))
}

# repetition function
monte.rep.rich.subs.moments <- function(parallel=TRUE,year.monte,vdgp){ #get equilibrium market shares and prices for each market and append them
  if(parallel==T)  foreach(i=1:n.reps, .combine = 'rbind')  %dorng%  monte.func.rich.subs.mlog.moments(i,year.monte,vdgp) 
  else foreach(i=1:n.reps, .combine = 'rbind')  %do%  monte.func.rich.subs.mlog.moments(i,year.monte,vdgp) 
}


### Function for computing CF for BLP data ###
#estimation of eta for BLP data
estapprox.BLPdata <- function(pre) { # using data on market shares and prices, costs and observable quality, 
  #estimate eta (and other parameters) under the CES approximation
  # CES market share #
  share.ces <- function(results) {
    exp(results$fitted.values)/sum(exp(results$fitted.values))
  }
  #
  mc <- pre$cit
  s <- pre$s
  p <- pre$p
  x <- pre$x
  f <- pre$f
  i <- pre$i
  # estimate the demand and quality coefficient ("CES" and "B")
  # 
    ols.results <- pre[,lm(log(s)~log(mc)+hpwt+air+mpd+space+factor(f)+xi)]
  #
  #
  cost.elas <- ols.results$coefficients[[2]]
  eta <- -cost.elas
  return(data.frame(pre,eta))
}



#CF function for BLP data
CF_BLPdata.simul <- function(vsim,cf.meth="EHA"){
  cat(paste0("Starting version: ",vsim,"\n \n"))
  cat(paste0("Initializing: ","\n"))
  setwd("~/Dropbox/BLP_REStat/Replication")
  #Initialization phase
  data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=3) #gen exog data with v=3 for initialization
  BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective=objective) 
  BLP.output.m <- BLP.output[m==j]
  BLP.meanownelas <<- BLP.output.m[, mean(own_elas_m)] # super assignment needed here because of scope for DGP.BLP
  BLP.priceadjust <<- BLP.output.m[, mean(p_m * (1-s_m))] 
  #  v=1 version --> wvp_m = 0 --> average own elasticity will be exactly same as v=3
  if(vsim<3){
    cat(paste0("Re-Initializing for version ",vsim,"\n"))
    set.seed(seedn)
    data.exog.mm <- DGP.BLP(year.chosen=year.monte,vdgp=vsim) #gen exog data
    BLP.output <- rich.subs.BLP.mlog(dat=data.exog.mm,verbose=FALSE,objective=objective) 
    BLP.output.m <- BLP.output[m==j]
    BLP.output.m[, own_elas.theory := -(alpha.base/ybar) * p_m * (1-s_m-wvp_m/s_m)]
    cat(paste0("Theoretical own elas match: ",BLP.output.m[, all.equal(-own_elas.theory,-own_elas_m)],"\n"))
    #note that for v = 1, wvp_m should be 0, and therefore line 74 still holds
    #alpha.base is outputed in DGP.BLP
  } 
  #
  #pre 
  Y <- sum(data.exog.mm$households$y.h) # total income of that economy
  cat("mean = ",BLP.output.m$lownelas.mean.hh[1],"sd = ",BLP.output.m$lownelas.sd.hh[1],file=paste0("Tables/CF_BLPdata/hh_elas_v",vsim,".txt"))
  #instead of estimating eta, we assume we already know it.
  pre <- BLP.output.m[,.(Y=Y,eta = -BLP.meanownelas,own_elas_m,super_elas_m,m,year,f,p=p_m,s=s_m,i,n,xi=xi.predicted_m,mc=mc_m,cit=mc_m,const,hpwt=hpwt_m,air=air_m,mpd=mpd_m,space=space_m)]
   #now add tariff change, xvars, mc and xi to original data
  n.models <- nrow(data.exog.mm$models)
  data.exog.mm$models <- data.exog.mm$models[,.(year,f,i,n
                                                ,m=1:n.models
                                                ,x0=const
                                                ,x1=hpwt
                                                ,x2=air
                                                ,x3=mpd
                                                ,x4=space
                                                ,xi=pre$xi
                                                ,mc=pre$mc
                                                ,tau=1
                                                ,tariff.change=tar.increase)]
  setnames(data.exog.mm$households, old=c("b.h.const","b.h.hpwt","b.h.air","b.h.mpd","b.h.space"),new=c("b0.h","b1.h","b2.h","b3.h","b4.h"))
  data.exog.mm$households <- cbind(data.exog.mm$households,n=market)
  #post
  setwd("~/Dropbox/BLP_REStat/Replication/Rfiles")
  cat(paste0("Computing Post ","\n"))
  post <- eqbm.mlog(dat=data.exog.mm,market=market,v=vsim,tar.change=1)
  setwd("~/Dropbox/BLP_REStat/Replication")
  #
  #PT
  pt.model <- (post$p[post$i!=post$n]-pre$p[pre$i!=pre$n])/(post$cit[post$i!=post$n]- pre$cit[pre$i!=pre$n])
  pte.model <- pt.model * (pre$cit[pre$i!=pre$n]/pre$p[pre$i!=pre$n])
  oce.model <- pre$own_elas_m[pre$i!=pre$n] * pte.model #own cost elasticity model level: own price elasticity times pte.model 
  pt <- mean(pt.model[pt.model!= -Inf & pt.model!= Inf]) # pass thru rate (IO)
  pte <- mean(pte.model[pte.model!= -Inf & pte.model!= Inf]) #  pass thru elasticity  a la Int'l Macro (IM)
  oce <- mean(oce.model[oce.model!= -Inf & oce.model!= Inf]) #  own cost elasticity 
  s.f <- pre$s
  s.f[pre$i==pre$n] <- 0 #among foreign models
  s.f[post$cit-pre$cit==0] <-0 # among foreign models that have some variation in tariffs.
  r.f <- rank(-s.f)
  pte1 <- ((post$p[r.f==1]-pre$p[r.f==1])/(post$cit[r.f==1]- pre$cit[r.f==1]))*(pre$cit[r.f==1]/pre$p[r.f==1]) #  pass thru elasticity  a la Int'l Macro (IM) for top foreign firm NOT A MEAN
  #
  pre <- within(pre,{ptm = pte}) # assume we know the pass-through elasticity
  #
  #VALUE share comparison
  #rescale market shares to be relative to inside good.
  pre <- within(pre,{s.true = s}) # quantity share of relative to all expenditures
  post <- within(post,{s.true = s}) # quantity share of relative to all expenditures
  pre <- within(pre,{s = p*s/sum(p*s)}) # value share of insiders (total quant cancels in s.true)
  post <- within(post,{s = p*s/sum(p*s)}) # value share of insiders
  #do CES approximation on value shares
  post.with.ces <- as.data.table(CF.approx(pre,post,meth=cf.meth)) #add ces approximation results
  #
  #Now proceEd to comparisons:
  setDT(pre)
  setDT(post)
  #
  #true vs approx changes
  chg.true <- round(100*(post.with.ces[i==n,sum(s)]-pre[i==n,sum(s)]),2)
  chg.ces  <- round(100*(post.with.ces[i==n,sum(s.ces)]-pre[i==n,sum(s)]),2)
  bias.value <- chg.ces-chg.true
  #
  #
  #QTY share comparison
  #reinstate CF.approx based on quantity shares on the inside good
  pre <- pre[, s := s.true/sum(s.true)] # qty share of insiders
  post <- post[, s := s.true/sum(s.true)] # qty share of insiders
  post.with.ces <- as.data.table(CF.approx(pre,post,meth=cf.meth)) #add ces approximation results
  #
  #true vs approx changes
  chg.true.qty <- round(100*(post.with.ces[i==n,sum(s)]-pre[i==n,sum(s)]),2)
  chg.ces.qty  <- round(100*(post.with.ces[i==n,sum(s.ces)]-pre[i==n,sum(s)]),2)
  bias.qty <- chg.ces.qty-chg.true.qty
  #below are diagnostics for understanding the fact that v2 outcome falls compared to v1
  chg.OG.qty <- post[,1-sum(s.true)]-pre[,1-sum(s.true)]
  saveRDS(pre,paste("Tables/CF_BLPdata/pre_setting",vsim,".rds",sep=""))
  saveRDS(post,paste("Tables/CF_BLPdata/post_setting",vsim,".rds",sep=""))
  #
  #Table of results
  CF.BLP.output <- data.table(v=vsim,eta=pre$eta[1],oce=round(oce, digits=2),meanownelas = round(mean(pre$own_elas_m),digits=2)
                              ,meansuperelas = round(mean(pre$super_elas_m),digits=2),pt=round(pt,digits = 2),pte=round(pte, digits=2),pte1=round(pte1,digits=2)
                              ,chg_true_qty=chg.true.qty,chg_ces_qty=chg.ces.qty,bias_qty=bias.qty
                              ,chg_true=chg.true,chg_ces=chg.ces,bias_value=bias.value,chg.OG.qty=chg.OG.qty)
  return(CF.BLP.output)
}

